/***************************

This function loads or sets the sfh
We read from the file in two columns (time in yr
and SFR in M/yr) and then we make a comulative
mass using this equation

Sum SFR_i*(t_{i+1}-t_i)


The first entry in SFR is put to zero.
Timesteps have to be sorted.

***************************/

//read in globals
extern float galaxy_sfr;
extern string galaxy_sfh;
extern float maxtime;


void setsfh(double **sfh_time, double **sfh_cumul,int *npntsfh){


  //check if continuous SFR or SFH
  if(galaxy_sfr > 0){

    //continuous SFR: we are building a constant SFH that will be used later
    //20 is arbitrary

    //reallocate malloc
    (*npntsfh)=20;
    *sfh_time = (double*) realloc(*sfh_time,20*sizeof(double));
    *sfh_cumul = (double*) realloc(*sfh_cumul,20*sizeof(double));
    for(int i=0;i<20;i++){
      (*sfh_time)[i]=(double) maxtime*(i/19.);
      (*sfh_cumul)[i]=(double) galaxy_sfr*(*sfh_time)[i];
    }
    
  } else {
    
    //use SFH from user
    ifstream infile;
    string line;
    
    infile.open(galaxy_sfh.data(),ifstream::in);
    if (infile.is_open())
      {
	//get number of lines 
	while (line.compare("END") != 0)
	  {
	    getline(infile,line);
	    (*npntsfh)++;
	  }  
	(*npntsfh)=(*npntsfh)-2;
	infile.close();
      } else {
      cout<<"SETSFH: I cannot find the file "<<galaxy_sfh<<endl;
      exit(2);
    }
    
  
    //reallocate malloc
    *sfh_time = (double*) realloc(*sfh_time,(*npntsfh)*sizeof(double));
    *sfh_cumul = (double*) realloc(*sfh_cumul,(*npntsfh)*sizeof(double));
    
    float buf,buf1;
    char bufc[20];
    FILE * pFile;
    pFile = fopen (galaxy_sfh.data(),"r");
    fscanf(pFile,"%s",bufc);
    for(int i=0;i<(*npntsfh);i++){
      fscanf(pFile,"%e %e",&buf,&buf1);
      if(i==0){
	//set first timestep
	(*sfh_time)[i]=(double) buf;
	(*sfh_cumul)[i]=(double) 0;
      } else {
	//set other timesteps
	(*sfh_time)[i]=(double) buf;
	(*sfh_cumul)[i]=(double) (buf1*(buf-(*sfh_time)[i-1]));
      }
    }
    fclose(pFile);
    
    //make the cumulative
    for(int i=(*npntsfh)-1;i>=0;i--){
      for(int j=i-1;j>=0;j--){
	(*sfh_cumul)[i]+=(*sfh_cumul)[j];
      }
    }
  }
  return;
}
